

####################################################################################
#required to build code OUTSIDE R
cd /home1/31/jc165798/SCRIPTS/sdmcode/misc/ben.phillips/individual.genetic.models/

R CMD SHLIB Functions.c


####################################################################################
#START R
#the function
dyn.load("/home1/31/jc165798/SCRIPTS/sdmcode/misc/ben.phillips/individual.genetic.models/Functions.so")
system.time( .Call('mother',R_n=10000,R_spX=10,R_spY=10,R_lambda=5.0,R_K=1000,R_disp_cost=0.1,R_reports_rate=20,R_num_gens=101,rho=environment(), R_m_rate=0.0001))

# R_n - starting population size
# R_spX & R_spY - starting habitat matrix size
# R_lambda - basic number of offspring
# R_K - carrying capacity
# R_disp_cost - cost of dispersing .. ranges between 0 & 1
# R_reports_rate - rate at which outputs of density & individuals are sent back to R
# R_num_gens - the number of generations past the burn in period... plus also the increase spX to be used
# rho - the R working environment






#### Postprocessing functions ####

#calculates linkage disequilibrium scaled against Dmax
# assumes data are genotypes in a matrix [individual][D1a, D1b,...Nna, Nnb]
# assumes diploid, bi-allelic system
# would pass either selected or neutral arrays to mtrx (not both!)
# returns mean R = D/Dmax across all locus pairs
ld<-function(mtrx, nloci, scale=TRUE){
	mtrx<-subset(mtrx, mtrx[,1]!=-9) #remove -9s
	gtype.pairs<-matrix(combn(1:nloci, 2), nrow=2) #find pairwise combinations of loci
	a.pairs<-matrix(seq(1, 2*nloci, by=2)[gtype.pairs], nrow=2) #a chromosome column pairs
	b.pairs<-matrix(seq(2, 2*nloci, by=2)[gtype.pairs], nrow=2) #b chromosome column pairs
	R<-c() #vector to take scaled D values
	for (ii in 1:length(gtype.pairs[1,])){ #step through locus pairs
		L1<-factor(mtrx[,a.pairs[1,ii]], levels=0:1)
		L2<-factor(mtrx[,a.pairs[2,ii]], levels=0:1)
		gam.fqs<-table(L1, L2) #creates table of P11, P12, etc.
		L1<-factor(mtrx[,b.pairs[1,ii]], levels=0:1)
		L2<-factor(mtrx[,b.pairs[2,ii]], levels=0:1)
		gam.fqs<-gam.fqs+table(L1, L2) #creates table of P11, P12, etc.
		gam.fqs<-gam.fqs/sum(gam.fqs) #as relative frequencies
		p1p2<-apply(gam.fqs, 1, sum) #vectors of allele frequencies
		q1q2<-apply(gam.fqs, 2, sum)
		if (sum(c(p1p2, q1q2)>0.9)>=1) next #avoid nonsense
		D<-gam.fqs[1,1]-p1p2[1]*q1q2[1]
		if (D>0 && scale==T) { #scale
			Dmax<-min(p1p2[1]*q1q2[2], p1p2[2]*q1q2[1])
			D<-D/Dmax	
		}
		if (D<0 && scale==T) { #scale
			Dmax<-max(-p1p2[1]*q1q2[1], -p1p2[2]*q1q2[2])
			D<-D/Dmax	
		}
		R<-c(R, D)
	}
	R<-mean(R, na.rm=T)
	return(R)		
}


#runs a transect through population from L to R and calculates LD along trnasect with a moving window
#	of scale (grid cells) equal to neigh.size
tsect<-function(mtrx, neigh.size, spY, nloci){
	mtrx<-subset(mtrx, mtrx[,1]!=-9) #remove -9s
	y.min<-1
	y.max<-neigh.size
	if (y.max>spY) return("Error - neigh.size cannot exceed spY")
	mtrx<-subset(mtrx, mtrx[,2]>=y.min & mtrx[,2]<=y.max) #trim mtrx by y
	x.min<-5
	max.occ.x<-max(mtrx[,1])
	tsect.sites<-seq(max.occ.x, x.min, -neigh.size)
	mtrx<-subset(mtrx, mtrx[,1]>=min(tsect.sites)) # trim by x
	tsect.sites<-rev(rep(tsect.sites, each=neigh.size)) # renumber x to reflect tsect.sites
	mtrx[,1]<-tsect.sites[mtrx[,1]-min(tsect.sites)+1]
	cols<-length(mtrx[1,-(1:2)])
	mtrx<-split(mtrx[,-(1:2)], mtrx[,1])
	mtrx<-lapply(mtrx, matrix, ncol=cols)
	out<-lapply(mtrx, ld, nloci=nloci)
	out 	
}




####### Test functions (redundant) #########

#function to transform matrix of alleles to genotype dataframe for LD
tf.data<-function(mtrx){
	cnames<-c()
	for (ii in 1:(length(mtrx[1,])/2)){
		vname<-paste("loc", ii, sep="")
		cnames<-c(cnames, vname)
		assign(vname, genotype(mtrx[,1+(ii-1)*2], mtrx[,2+(ii-1)*2]))
		if (ii==1) out<-data.frame(vname=get(vname))
		if (ii>1) out<-cbind(out, vname=get(vname))
	}
	colnames(out)<-cnames	
	out
}


cond.LD<-function(dframe){
	gtype.prs<-matrix(combn(1:length(dframe[1,]), 2), nrow=2)
	out<-c()
	for (ii in 1:length(gtype.prs[1,])){
		temp<-dframe[, gtype.prs[,ii]]
		if (length(levels(temp[,1]))==1 | length(levels(temp[,2]))==1){
			out<-c(out, NA)	
		}	
		else {
			out<-c(out, LD(temp)[[3]][1,2])	
		}
	}
	out	
}

tsect2<-function(mtrx, neigh.size, spY, nloci){
	require(genetics)
	mtrx<-subset(mtrx, mtrx[,1]!=-9) #remove -9s
	y.min<-1
	y.max<-neigh.size
	if (y.max>spY) return("Error - neigh.size cannot exceed spY")
	mtrx<-subset(mtrx, mtrx[,2]>=y.min & mtrx[,2]<=y.max) #trim mtrx by y
	x.min<-5
	max.occ.x<-max(mtrx[,1])
	tsect.sites<-seq(max.occ.x, x.min, -neigh.size)
	mtrx<-subset(mtrx, mtrx[,1]>=min(tsect.sites)) # trim by x
	tsect.sites<-rev(rep(tsect.sites, each=neigh.size)) # renumber x to reflect tsect.sites
	mtrx[,1]<-tsect.sites[mtrx[,1]-min(tsect.sites)+1]
	cols<-length(mtrx[1,-(1:2)])
	mtrx<-split(mtrx[,-(1:2)], mtrx[,1])
	mtrx<-lapply(mtrx, matrix, ncol=cols)
	mtrx<-lapply(mtrx, tf.data)
	out<-lapply(mtrx, cond.LD)
	out 	
}






